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Abstract 

We  describe  a  second-order  finite- volume  formulation  for  solving  incompressible  flow 
problems  on  Cartesian  grids  using  a  colocated  arrangement  of  variables.  A  projec¬ 
tion  method  is  used,  in  which  a  provisional  cell-centered  velocity  is  obtained  by 
integrating  the  effects  of  advection  and  viscous  forces  over  the  numerical  time  step. 
Interface  velocities  are  then  interpolated  from  the  provisional  cell-centered  velocities 
and  corrected  using  the  solution  of  a  Poission  equation  for  the  pressure  distribution. 
The  cell-center  velocities  are  then  interpolated  from  the  corrected  interface  veloci¬ 
ties.  This  pair  of  interpolation  steps  adds  numerical  diffusion  that  mimics  a  normal 
viscous  stress,  and  the  form  of  this  damping  term  is  similar  to  that  which  comes 
from  stabilizing  the  Poisson  equation  for  pressure  by  adding  a  term  that  is  propor¬ 
tional  to  the  fourth-derivative  of  pressure.  The  method  requires  no  extra  damping 
terms  in  order  to  maintain  stability. 
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Manuscript  approved  December  10, 2004. 


1  Introduction 


Many  engineering  applications  involve  incompressible  flows,  from  pipe-flow 
problems  to  microfluidics  applications.  A  wide  range  of  finite-element,  finite- 
difference,  and  finite-volume  methods  for  solving  these  problems  numerically 
have  been  developed  over  the  years.  [1-9].  Many  of  the  more  popular  finite- 
volume  schemes  are  clasified  as  pressure-correction  methods.  In  these  meth¬ 
ods,  a  provisional  velocity  field  for  the  new  time  level  or  iteration  in  a  steady 
state  calculation  is  calculated  that  includes  the  effects  of  viscous  and  advective 
transport  but  not  the  full  effect  of  the  pressure  gradient.  A  Poisson  equation 
for  the  new  pressure  distribution  (or  for  a  correction  to  the  pressure  distri¬ 
bution)  is  then  solved,  and  this  solution  enforces  mass  conservation  once  the 
velocity  distibution  is  updated  using  the  new  pressures.  The  well-known  SIM¬ 
PLER],  SIMPLEC[4],  QUICK[5],  and  PISO[6]  methods  fall  into  this  category. 
Fractional  step  methods  can  be  grouped  with  these  methods,  but  distinctions 
between  various  fractional  step  methods  are  often  made  depending  on  if  and 
how  the  old  pressure  field  is  used  to  determine  the  provisional  velocity  field. 
The  unifying  theme  in  these  projection  methods  is  the  solution  of  a  Poisson 
equation  for  the  pressure  or  a  pressure-like  variable  that  provides  a  correction 
to  the  velocity  distribution  that  enforces  conservation  of  mass. 

The  application  of  SIMPLE-type  methods  was  first  successful  on  staggered 
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meshes  in  which  the  control  volumes  for  the  pressure  and  various  velocity 
components  were  offset.  This  placed  velocity  components  on  the  faces  of  the 
pressure  control  volume  such  that  central-differencing  of  the  pressure  Poisson 
equation  led  to  compact,  consistent  stencils.  Staggering  the  meshes,  however, 
leads  to  added  complexity  in  the  grid,  particularly  when  modeling  flows  around 
complex  geometries. 

If  pressure  and  velocity  variables  are  colocated,  one  set  of  control  volumes 
can  be  used  for  all  variables.  Now,  however,  a  different  problem  emerges.  The 
pressure  Poission  equation  used  in  projection  methods  is  derived  by  taking 
the  divergence  of  the  momentum  equation.  Therefore  the  Laplacian  of  pres¬ 
sure,  V2p,  in  the  pressure  equation  is  the  divergence  of  the  pressure  gradient 
term  from  the  momentum  equation,  V  ■  (Vp).  Once  the  discretizations  for  the 
divergence  and  gradient  are  chosen,  only  one  discretization  for  the  Laplacian 
guarantees  that  V2p  =  V  •  (Vp).  Choosing  this  discretization  ensures  that 
the  resulting  pressure  equation  accurately  imposes  V  •  v  =  0.  If  standard 
central-difference  stencils  are  used  to  discretize  the  divergence  and  gradient 
operators,  the  resulting  stencil  for  V  •  V  =  V2  produces  odd-even  uncoupling. 
That  is,  V2(f>  at  cell  (i,j)  uses  the  values  of  0  two  cells  away  from  (i,j)  (such 
as  ( i  +  2 ,j)  and  (i,j  —  2))  but  not  the  values  of  immediate  neighbors  (such 
as  (i  +  l,j)  or  (i,  j  + 1)).  Ignoring  the  requirement  that  the  discrete  operators 
satisfy  V2^>  =  V  •  (V0)  and  using  compact  stencils  for  all  three  operators  can 
produce  conservation  errors  and  instabilities.  Standard  ways  to  overcome  this 
instability  include  adding  a  damping  term  proportional  to  V4p  to  the  Pois- 
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son  equation  for  the  pressure  or  by  correcting  the  interface  velocities  used  to 
calculate  V  -v  [1,2,9], 

This  paper  describes  a  compact  finite- volume  method  for  solving  incompress¬ 
ible  flow  problems  on  a  colocated  grid.  Rather  than  use  the  pressure  distribu¬ 
tion  to  correct  the  cell-centered  values  directly,  the  normal  velocity  component 
at  each  cell  interface  (interpolated  from  the  provisional  cell-centered  values) 
is  corrected,  and  the  cell-centered  values  are  interpolated  from  the  interface 
values.  This  process  introduces  numerical  diffusion  that  has  the  same  form  as 
a  normal  viscous  stress  and  is  similar  to  the  popular  V4p  damping  approach 
for  stabilizing  the  pressure  equation.  Therefore,  no  extra  damping  terms  are 
included  in  the  Poisson  equation  for  the  pressure  in  order  to  ensure  stability. 


2  Governing  Equations  and  Basic  Discretization 


We  begin  with  the  incompressible  Navier  Stokes  equations, 


dpv 

dt 


+  V  •  (pvv)  +  Vp  —  /iV2v  =  0, 
V  •  v  =  0, 


(1) 

(2) 


in  which  p  is  the  fluid  density  and  p  is  the  viscosity.  The  constraint  V  •  v  =  0 
on  the  velocity  v  only  requires  that  the  substantial  derivative  of  p  be  zero: 


Dp  dp 

m  =  m+v'vp  =  0- 


(3) 
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We  make  the  additional  assumption  that  p  and  p  are  homogeneous  throughout 
the  computational  domain.  The  component  of  the  convective  term  V  •  (pvv) 
appearing  in  the  momentum  equation  for  direction  i  is  V  •  ( pupv ). 

Our  focus  is  the  spatial  discretization  of  the  flow  variables  and  differential 
operators,  so  we  will  consider  steady-state  calculations  obtained  by  marching  a 
simple  time-accurate  scheme  to  equilibrium.  This  basic  approach  allows  for  two 
paths  of  development  in  the  future:  improved  steady-state  convergence  rates  by 
implementing  a  local-time-stepping  algorithm,  and  improved  time  resolution 
for  unsteady  calculations  by  increasing  the  order  of  the  time  discretization. 
A  forward  Euler  time  discretization  between  times  tn  and  tn+i  separated  by 
time  step  At  yields 


pvn+1  -  pvn 


Ai  (  —  V  •  (pvv)n  -  Vpn+1  +  pV2v") 


(4) 


In  Eq.  (4),  we  have  treated  the  pressure  term  implicitly  and  the  advective  and 
viscous  terms  explicitly.  A  provisional  velocity  v  can  be  calculated  using 


pv  =  pvn  +  At  (  -  V  •  (pw)n  +  pV2vn )  (5) 

Now  replace  the  time-level  n  terms  in  Eq.  (4)  using  v  from  Eq.  (5)  and  take 
the  divergence  of  the  result.  Requiring  V  •  v  =  0  then  gives  the  differential 
form  of  the  pressure  equation: 

vy+i  =  -W  •  (/#)■  (0 
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Finally,  vn+1  is  calculated  using 


(pv)n+1  =pv-  AtVpn+l.  (7) 

From  the  derivation  of  Eq.  (6),  the  Laplacian  term  arises  from  composition  of 
the  divergence  operator  with  the  gradient  operator.  The  methods  that  com¬ 
pensate  for  the  odd-even  instability  either  add  a  term  proportional  to  V4p 
into  Eq.  (6),  or  correct  the  interpolated  interface  velocity  using  centered  and 
interpolated  pressure  gradients  [1].  They  therefore  solve  an  altered  form  of 
Eq.  (6)  to  find  the  pressure  distribution  and  then  use  Eq.  (7)  to  update  the 
velocity  values. 


3  Current  Method 

The  algortihm  is  presented  for  two-dimensional  flows,  but  the  extention  to 
three-dimensions  is  straightforward.  Three-dimentional  simulations  of  microchan- 
nel  flows  using  this  algorithm  have  been  presented  elsewhere  [12, 13].  The  spa¬ 
tial  discretization  is  done  on  a  Cartesian  grid,  as  illustrated  in  Figure  1.  Integer 
subscripts  indicate  cell  centered  locations,  and  fractional  subscripts  indicate 
interface  locations.  For  example,  pitj  is  the  pressure  at  point  and  ui+ 1/2  is 

the  x-velocity  at  the  interface  between  cells  (i,  j)  and  (i  +  1,  j).  In  the  current 
work,  the  variables  are  assumed  to  vary  linearly  within  each  cell,  so  the  point 
value  at  the  cell  center  and  the  cell-average  value  are  identical.  The  cell  acts  as 
the  control  volume  for  all  variables.  The  viscous  term  in  Eq.  (5)  is  discretized 
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using  a  flux  formulation  equivalent  to  second-order  centered  finite-differencing 
on  cell-centered  values.  The  effect  of  the  advection  term  in  Eq.  (5)  is  calculated 
using  flux-corrected  transport  (FCT)  as  implemented  in  the  LCPFCT  library 
[10, 11].  LCPFCT  solves  generalized  conservation  laws,  but  Eq.  (6)  specifies 
the  pressure  field  using  an  algebraic  constraint  rather  than  the  integration  of 
a  time  rate-of-change.  LCPFCT  is  used  to  calculate  the  contribution  of  the 
advection  term,  V  •  (pw)n,  to  v  so  that  Eq.  (6)  can  be  solved  for  the  updated 
pressure  field. 

Rather  than  write  Eq.  (6)  using  the  cell-centered  provisional  velocities,  we 
calculate  provisional  interface  normal  velocities  using  a  linear  interpolation. 
For  unifiorm  Ax  and  Ay,  this  gives 

P^i+l/2  —  2  T  Pui+l,j^j 

PVj+1/2  =  £  (pSij  +  fiUiJ+lj  ■  (8) 

Then  we  use  these  interpolated  normal  velocities  to  calculate  V  •  (pv): 

,  Flux  x  area 

V  ■  </,V)«  =  cell  volume 

=  AxAij  [aj/(|»5.+1/2  -  pUi-1/2)  +  Ax(pvj+i/2  -  pv,-,/2)  j ,  (9) 

in  which  the  sum  is  taken  over  the  interfaces  of  the  cell. 

We  then  solve  Eq.  (6)  using  a  flux-based  approach  that  is  equvalent  on  a 
regularly-spaced  grid  to  using  a  compact  central-difference  stencil  for  the 
Laplacian  term,  and  we  do  not  add  a  stabilizing  term  to  the  equation.  Rather 
than  update  the  cell-centered  values  directly  using  the  new  pressure  distribu- 
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tion,  we  correct  the  interface  normal  velocities  using 

(p«)Sl/2  =  (p«)i+ 1/2  -  At"+h&x  Pb3  > 

(H?£v2  =  (/w)i+i/2  -  (10) 

The  cell  centered  values  are  then  linearly  interpolated  from  the  interface  val¬ 
ues.  For  constant  Ax  and  Ay,  this  reduces  to 


Therefore,  the  solution  procedure  is 

(1)  Calculate  cell-centered  pv  using  Eq.  (5). 

(2)  Calculate  Ui+1/2 ,  #7+1/2  by  interpolating  the  cell-centered  values. 

(3)  Use  Ui+ 1/2,  £>7+1/2  to  calculate  V  •  (pv)  as  in  Eq.  (9). 

(4)  Solve  Eq.  (6)  for  the  cell-centered  pressure  values,  using  a  compact  stencil 
for  V2. 

(5)  Update  the  interface  normal  velocity  values  using  Eq.  (10). 

(6)  Update  cell-centered  velocities  using  Eq.  (11). 

The  difference  between  this  method  and  others  in  the  literature  is  the  use  of 
two  interpolation  steps:  first  to  calculate  provisional  interface  normal  veloc¬ 
ities  from  cell-centered  values,  and  the  second  to  calculate  the  cell-centered 
values  from  the  corrected  interface  velocities.  The  second  interpolation  step  is 
needed  because  the  interface  velocities,  and  not  the  cell-centered  velocities,  are 
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corrected  using  the  pressure  distribution.  This  interpolation  produces  a  stabi¬ 
lizing  effect,  so  that  there  is  no  need  to  add  stabilizing  terms  to  the  pressure 
equation. 


4  Numerical  Diffusion 


We  now  show  that  interpolating  the  normal  velocity  onto  interfaces  and  then 
back  to  cell  centers  produces  a  numerical  damping  term  that  has  the  same 
form  as  a  normal  viscous  diffusion.  Consider  a  function  u(x,y,t )  for  which 
values  on  cell  centers  Ui,U{+i  and  interfaces  rq+i/2  on  an  equally  spaced  grid 
are  known.  The  interpolated  value  ui+ x/2  is  given  by 


_  Hi  “f"  Hi- {_! 

ui+ 1/2  =  - 7) - 


(12) 


Then  interpolating  cell-centered  values  to  interfaces,  and  then  back  to  cell 
centers,  gives  (using  a  double  overline  to  denote  the  two  interpolation  steps) 


Ui 


Hi- 1/2  +  Uj+1/2 
2  . ~ 


1 

2 


+  ui )  +  2^Ui 


(13) 

(14) 

(15) 


with  52u/5x 2  indicating  a  centered  finite-difference  approximation  for  the  sec¬ 
ond  derivative.  This  pair  of  interpolations  occurs  once  per  timestep,  so  we  can 
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J 


view  this  diffusion  as  an  unsteady  term  by  rearranging  Eq.  (15)  as 


pui  —  pui  p/S.x2  ( 82u 

A t  4A t  l  8x2 


(16) 


The  discretized  viscous  term  in  the  ^-momentum  equation  is 


.  2  .  52u  82u  52u 


(17) 


so  the  interpolation  adds  numerical  diffusion  to  the  normal  component  only: 


((. iV2u)  +  numerical  diffusion 
,  .82u  82u  82u 

=  0‘  +  iv)^  +  /‘^  +  ^- 


(18) 


Defining  a  viscous  cfl  number,  cfl„  as 


cflv 


A  tp 
A  x2p' 


(19) 


the  numerical  diffusion  coefficient  pn  can  be  written  as 


V'Tl 


4cflv  * 


(20) 


This  is  an  approximate  analysis;  the  solution  to  Eq.  (6)  is  used  to  update  the 
interface  velocities  bewteen  the  two  interpolation  steps,  whereas  the  deriva¬ 
tion  of  Eq.  (20)  lumps  the  two  interpolation  steps  together  and  neglects  the 
changes  that  the  pressure  update  imposes.  The  pressure  update  imposes  a 
constraint  rather  than  calculating  an  evolution  step  based  on  a  source  term, 
so  the  pressure  calculation  does  not  necessarily  preserve  the  diffusion  that  re¬ 
sults  from  the  initial  interpolation.  Equation  (20),  however,  does  show  that 
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increasing  the  timestep  relative  to  the  viscous  timescale  will  reduce  the  effect 
of  this  numerical  diffusion. 

The  impact  of  the  numerical  diffusion  on  solution  acuracy  can  be  assessed  in 
two  ways:  the  magnitude  of  /in  relative  to  fi,  which  is  given  by  the  viscous  cfl 
number,  or  the  magnitude  of  ^nd2u/dx2  compared  to  the  (physical)  viscous, 
pressure,  and  advection  terms  for  each  of  the  coordinate  directions.  If  fj,n  <C  //, 
then  the  numerical  diffusion  will  not  degrade  the  accuracy  of  the  calculation. 
The  inverse  is  not  necessarily  true:  a  large  value  of  /in  compared  to  //  does  not 
guarantee  that  the  simulation  results  are  tainted  by  excess  diffusion.  In  many 
flows,  the  normal  viscous  stress  is  of  secondary  importance  compared  to  the 
effects  of  shear  stresses,  pressure,  and  advection.  Then  this  diffusion  would 
not  be  problematic  using  an  explicit  solver  with  a  modest  cfl  constraint  for 
such  cases.  If  the  global  time-step  constraint  imposed  by  advection  is  stricter 
than  that  for  viscous  transport,  it  could  be  advantageous  to  use  a  local  time- 
step  algorithm.  This  would  allow  each  cell  to  march  toward  steady-state  based 
on  the  local  value  of  the  velocity  and  would  increase  the  local  time  step  and 
reduce  the  effect  of  this  numerical  diffusion.  On  the  other  hand,  if  the  time- 
step  constraint  for  advection  is  less  stringent  than  for  the  viscous  transport,  an 
implicit  discretization  for  the  viscous  terms  would  ensure  that  this  numerical 
diffusion  is  insignificant. 

We  can  also  compare  the  effects  of  the  interpolation  to  the  effects  of  adding 
explicit  damping  terms  to  the  pressure  Poisson  equation.  Denoting  ui+ 1/2  as 
the  true  provisional  interface  velocity,  the  interpolated  provisional  velocity 
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vy+1  =  -J-v  •  (/»v) 

A  t  v  ' 

pAx 2  d3wn+1  pAj/2  d3un+1 
8At  5rr3  8A t  dyz 


Ax2  d4p 
+~8~a^  + 


Ay 2  94p 
8  dyA 


+  0(Ax4)  +  0(Ay4). 


(25) 


Therefore,  using  the  finite-difference  approach  to  calculate  the  divergence  in 
Eq.  (6)  using  interpolated  interface  velocities  is  equivalent  to  adding  two  sets  of 
terms  to  Eq.  (6).  The  first  set  is  a  third-derivative  of  each  velocity  component 
(as  one  might  expect  since  we  took  the  divergence  of  a  second-order  error  term 
in  the  provisional  velocity),  and  the  second  is  a  fourth-derivative  of  p  similar 
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to  the  damping  term  often  added  to  enhance  stability.  Therefore,  the  process 
of  interpolating  the  cell-centered  velocities  to  the  interfaces  mimics  in  part  the 
effect  of  including  a  fourth-derivative  of  pressure  as  a  damping  term  in  the 
pressure  equation. 

Finally,  an  alternate  way  of  viewing  the  interpolation  is  to  consider  each  in¬ 
terpolation  step  as  a  mapping:  one  mapping,  B,  from  cell- aver  aged  values 
to  point  values  at  the  interfaces,  and  a  second,  D,  from  point  values  to  cell 
averages.  If  these  mappings  were  perfectly  matched,  then  DB  =  I,  that  is, 
mapping  to  the  interfaces  and  then  back  again  would  produce  the  original 
data  set.  As  a  minimal  requirement,  the  mappings  must  be  conservative:  if 
4>'  =  DB(f) ,  then  ||0'||  =  |j^||  for  conserved  quantity  (j>  and  the  1-norm  ||  •  ||  that 
measures  the  amount  of  (j>  in  the  domain.  The  current  method  satisfies  this  re¬ 
quirement:  the  interpolation  steps  are  diffusive,  but  they  are  also  conservative. 
Additional  methods  can  be  pursued  by  replacing  the  linear  interpolation  used 
in  the  current  work  with  other  mappings  between  cell-centered  and  interface 
values. 


5  Numerical  Results  and  Discussion 

The  current  algorithm  has  been  implemented  in  the  code  TINY3D,  and  two- 
dimensional  examples  are  shown  here  to  demonstrate  the  method.  Kaplan  et 
al.  [12, 13]  present  three-dimensional  calculations  for  a  channel  micromixer 
including  a  herringbone-patterned  surface  geometry. 
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Figure  2  shows  the  flowfield  in  a  two-dimensional  channel  100  pm  high  and 
400  pm  long.  The  flow  on  the  left  enters  at  a  uniform  speed  of  1  cm/s,  and 
p  =  0.014  gm/(cm-s).  We  enforce  no-slip  conditions  on  the  top  and  bottom 
walls,  and  the  outflow  (right)  boundary  is  held  at  p  =  1  atm.  Due  to  the  uni¬ 
form  inflow,  a  short  development  region  is  seen  during  which  the  equilibrium 
parabolic  profile  is  formed.  After  this  developing  region,  the  pressure  drop  be¬ 
comes  linear.  As  is  seen  in  Table  1,  the  calculated  values  for  the  equilibirum 
pressure  drop  are  converging  to  the  theoretical  value  as  Ax  — >  0.  Figure  3 
illustrates  that  the  code  demonstrates  second-order  convergence  for  this  cal¬ 
culation.  The  simulation  using  40  x  160  cells  gives  dp/dx  =1678  dyne/cm3. 
This  is  within  0.1%  of  the  theoretical  value  of  1680  dyne/cm3.  For  this  case, 
the  damping  term  introduced  by  the  interpolation  is  of  no  consequence  since 
no  normal  velocity  gradients  exist  in  the  equilibrium  profile.  The  curve  in 
Fig.  3  is  the  parabola  that  crosses  the  y-axis  at  the  theoretical  value  of  1680 
dyne/cm3  and  passes  through  the  computational  data  point  at  Ax  =  2.5 pm. 
The  fact  that  this  same  parabola  passes  through  the  other  two  computational 
data  points  at  Ax  =  5 pm  and  Ax  =  10pm  indicates  that  the  error  is  propor¬ 
tional  to  Ax2  and  the  method  exhibits  Ax2  convergence. 

A  second  set  of  calculations  focused  on  the  two-dimensional  flow  over  a  backward¬ 
facing  step.  The  geometry,  x- velocity  distribution,  and  streamline  traces  (path¬ 
lines)  are  shown  in  Figure  4.  The  step  consists  of  a  1  mm  channel  opening  into 
a  1.94  mm  channel,  and  the  computational  domain  extends  10  mm  upstream 
(left)  of  the  step  and  20  mm  downstream.  Water  at  13  cm/s  enters  from  the 
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left  with  /j,  =  0.0089  cm2/s.  Following  Armaly  et  al.  in  their  seminal  work 
on  the  backward-facing  step  problem  [14],  we  calculate  a  Reynolds  number  of 
292  using  this  mean  velocity  as  the  characteristic  speed  and  twice  the  1  mm 
entrance  width  as  the  characteristic  length.  At  this  Reynold  number,  a  single 
vortex  forms  downstream  from  the  step.  The  size  of  this  vortex,  as  charac¬ 
terized  by  the  reattachment  point  on  the  bottom  wall,  provides  a  test  for 
validating  the  order  and  the  accuracy  of  the  code.  In  order  to  match  the  exact 
geometry  using  moderate  resolution,  variable  grid-spacing  is  used.  The  top  of 
the  step  corresponds  to  y  —  0.094  cm.  Cells  have  uniform  Ax  throughout  the 
grid,  Ay  =  Ax  for  those  cells  that  lie  above  y  =  0.094  cm,  and  Ay  —  0.94Aa; 
for  cells  that  lie  below  y  =  0.094  cm.  Therefore,  if  N  cells  span  the  entrance 
channel,  then  2N  cells  span  the  exit  channel,  and  the  dimension  of  the  two 
channels  match  the  1:1.94  ratio  from  the  experiment. 

This  geometry  and  set  of  flow  conditions  provides  a  severe  test  of  the  numerical 
diffusion  in  the  current  explicit,  time-accurate  implementation,  as  character¬ 
ized  by  Eqs.  (18)  and  (20).  The  time  step  for  the  calculation  is  limited  by  the 
advection,  so  cfl„  is  low,  and  there  is  significant  normal  viscous  stress  in  the 
region  just  downstream  from  the  step.  Figure  5  indicates  the  relative  magni¬ 
tudes  of  the  normal  and  transverse  x-velocity  derivatives  that  contribute  to 
the  viscous  term  for  the  x-momentum  equation.  The  ratio  shown  is 

_  \\d2u/dx2\\ 

\\d2ufdx2\\  -|-  \\d2u/dy2\\ 

For  much  of  the  computational  domain,  this  ratio  is  less  than  0.1,  indicating 
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that  the  effects  of  shear  are  at  least  nine  times  larger  than  the  normal  stress. 
There  are,  however,  some  regions  in  which  the  normal  viscous  stress  is  greater 
than  the  tangential  stress  (indicated  by  r  >  0.5),  most  noticeably  just  down¬ 
stream  of  the  step.  A  close-up  of  this  region  is  shown  in  Fig.  (6).  At  steady 
state,  the  net  source  term  for  each  cell  is  zero,  so  that  the  contributions  due 
to  viscous  flux,  advective  flux,  and  pressure  are  balanced.  Figure  (6)  shows 
the  individual  cells  for  the  medium-resolution  simulation.  The  value  of  r  for 
each  cell  indicates  the  relative  magnitude  of  the  (physical)  normal  stress  and 
shear  stress  acting  on  the  x- velocity  component  at  equilibrium.  High  values  of 
r  indicate  where  the  shear  stresses  on  the  x-velocity  do  not  dominate  normal 
stresses  for  this  same  velocity  component.  If  r  is  low  throughout  the  domain 
and  cfl„  is  high,  we  know  that  the  numerical  diffusion  introduced  by  the  pair 
of  interpolation  steps  will  be  insignificant.  In  this  region  with  high  values  of 
r,  the  y-velocity  is  dominant,  and  Figure  (6)  does  not  indicate  whether  the 
normal  stresses  on  the  x-velocity  component  play  a  significant  role  in  the  over¬ 
all  flow  solution.  Given  the  high  values  of  r  and  the  the  small  values  of  cfl„, 
the  calculation  must  be  compared  with  experimental  data  to  demonstrate  the 
impact  that  the  interpolation  steps  have  on  the  accuracy.  One  is  tempted  to 
compare  ixn\\d2u/dx2\\  to  the  pressure  and  advective  terms,  but  as  noted  in 
the  previous  section,  such  a  comparison  would  be  inconclusive.  The  test  of 
accuracy  for  this  calculation  must  be  a  study  of  the  convergence  behavior  of 
the  method  and  a  comparison  with  experimental  data. 

To  test  the  convergence  rate  and  accuracy,  a  series  of  calculations  on  various 
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grids  were  performed.  Table  2  gives  the  reattachment  lengths  for  several  grid 
sizes,  and  these  values  are  plotted  in  Figs.  7  and  8.  At  Re  =  292,  Armaly  et  al. 
reports  a  reattachment  length  normalized  by  the  height  of  the  step  of  6.55  [14], 
and  this  value  is  marked  in  Fig.  7.  The  agreement  between  the  computational 
results  and  the  experiment  is  very  good,  even  though  the  values  of  cfl„  were 
very  low  because  the  time  step  for  these  cases  is  limited  by  the  advective 
and  not  the  viscous  time  scale.  The  standard  cfl  number  for  advection  for  a 
characteristic  speed  w  is 


cfl  = 


wAt 

Ax 


(27) 


To  guarantee  positivity,  FCT  requires  cfl  <  0.5,  so  the  limit  for  the  time  step  is 
set  to  cfl  =  0.45,  which  results  in  a  more  stringent  constraint  than  limiting  the 
timestep  based  on  cfl,,.  The  resulting  values  of  cfl„  are  also  listed  in  Table  2. 


The  coarsest  grid  included  in  Fig.  (7)  used  only  five  cells  to  span  the  inlet  chan¬ 
nel  and  five  to  span  the  height  of  the  step.  At  this  resolution  a  single  vortex 
is  predicted  downstream  of  the  step,  but  the  computation  does  not  correctly 
capture  the  basic  shape  of  this  vortex.  The  other  four  calculations  correctly 
capture  the  shape  of  the  vortex  as  illustrated  in  Fig.  (4  ),  so  these  computations 
are  used  to  study  the  convergence  characteristics  of  the  algorithm.  The  results 
for  these  four  cases  are  shown  in  Fig.  (8).  Figure  (8)  also  includes  an  ideal  curve 
for  second-order  convergence  similar  to  the  theoretical  curve  in  Fig.  (3).  For 
the  channel-flow  calculations,  the  theoretical  value  of  the  steady-state  pressure 
gradient  is  known,  and  the  pressure  gradient  predicted  by  the  code  converged 
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to  the  theoretical  value  as  Ax  — >•  0.  A  theoretical  value  for  the  reattachment 


length  is  not  known,  and  the  experimental  value  cannot  serve  as  a  substitute 
for  this  value  due  to  experimental  limitations  and  uncertanties.  Therefore,  the 
curve  in  Fig.  (8)  indicating  ideal  second-order  convergence  behavior  was  gen¬ 
erated  using  the  data  points  for  the  two  smallest  grid  sizes.  Using  these  two 
data  points  and  extrapolating  to  Ax  =  0,  second-order  convergence  predicts  a 
reattachment  length  of  6.563.  If  we  define  E(Ax)  as  the  difference  between  the 
reattachment  length  calculated  using  grid  size  Ax  and  the  asymptotic  value  as 
Ax  — y  0,  then  a  second-order  method  will  give  E( Ax/2)  =  1/4  E{ Ax).  As  the 
graph  illustrates,  the  result  on  the  grid  with  Ax  =  0.1  mm  is  well  below  the 
theoretical  curve.  This  indicates  that  the  code  is  converging  at  a  rate  faster 
than  second-order:  i?(0.05mm)  ss  1/40  .E/0. 1mm).  Since  the  viscous  and  ad- 
vection  transport  routines,  interpolations,  and  pressure  distribution  and  gra¬ 
dient  calculations  are  second-order,  this  apparent  increase  in  the  convergence 
rate  is  due  to  the  reduced  diffusion  in  the  interpolation  of  cell-centered  to 
interface  quantities  as  n  decreases.  Not  only  are  the  operator  discretizations 
and  flux  calculations  converging  as  Ax2,  but  these  operators  also  are  working 
with  less-diffused  data. 


6  Summary 

A  method  for  interpolating  between  cell-centered  and  interface  values  leads 
to  a  straightforward  application  of  a  projection  method  for  incompressible 
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flows  on  colocated,  Cartesian  meshes.  The  interpolation  procedure  mimics,  in 
part,  a  standard  stabilizing  term  that  is  proportional  to  the  fourth-derivative 
of  pressure,  so  the  current  method  does  not  include  an  explicit  stabilizing 
term  in  the  pressure  equation.  Analysis  indicates  that  the  numerical  diffusion 
introduced  by  the  interpolation  steps  acts  as  a  normal  stress  and  increases 
as  the  ratio  of  the  numerical  time  step  to  the  viscous  time  scale  decreases. 
Comparisons  with  experimental  results  for  two-dimensional  backward-facing 
step  flow  demonstrate  that  accuracy  is  maintained  even  when  the  time  step  is 
limited  by  advection  and  is  much  greater  than  the  times  step  allowed  by  the 
viscous  transport. 
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Table  1 


Pressure  gradient  results  for  the  channel  flow  using  three  different  grids.  The  theo¬ 
retical  value  is  1680  dyne/cm3. 


Number  of  Cells 

dp/ dx 

%  error 

(height  x  length) 

dyne/cm3 

10  x  40 

1647 

1.96% 

20  x  80 

1672 

0.48% 

40  x  160 

1678 

0.12% 
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Table  2 


Normalized  reattachment  length  downstream  of  the  step  for  1:1.94  step  at  Re  = 
292,  for  various  grid  sizes.  Experimental  value  for  the  1:1.94  channel  is  6.55  [14]. 
Simulation  time  step  is  specified  by  the  requirement  that  cfl  =  0.45  for  the  advection 
routine.  The  error  is  relative  to  the  experimental  value. 


Cells 

spanning  inlet 

cflfl 

Ratio  of  reattachment 

length  to  step  height 

% 

error 

5 

0.012 

2.553 

61.1% 

10 

0.021 

5.620 

14.2% 

15 

0.029 

6.470 

1.3% 

20 

0.039 

6.539 

0.2% 

40 

0.080 

6.557 

0.1% 
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Fig.  1.  Grid  arrangement  for  the  current  method.  All  variables  are  stored  at  cell 
centers,  and  the  cell  acts  as  the  control  volume  for  all  variables. 
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2.  X-velocity  distribution  for  the  2D  channel  test  problem. 
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Fig.  3.  Steady-state  pressure  gradient  in  the  100pm  microchannel  as  a  function  of 
grid  size.  Open  circles  are  computational  results;  filled  circle  is  the  exact  solution. 
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Fig.  4.  x-velocity  contours  and  stream  traces  for  the  backward-facing  step  problem. 
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Fig.  5.  The  ratio  \\d2u/ dx2\\/ (\\d2u/ dx2\\  +  ||<32u/ch/2||),  which  indicates  the  relative 
importance  of  the  normal  viscous  stress  compared  to  the  tangential  viscous  stress 
at  steady  state.  Contours  generated  by  treating  the  cell-centered  values  as  point 
values. 


28 


(Uio)  A 


r:  o.i  0.5  0.9  5 


X  (cm) 


Fig.  6.  Detail  behind  the  step  of  the  ratio  \\d2u/ dx2\\/ {\\d2u!  dx2]^  +  \\d2u/dy2\\)  at 
steady  state,  plotted  as  cell  values. 
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Fig.  7.  Normalized  reattachment  length  as  a  function  of  grid  size.  Circles  are  com¬ 
putational  results;  dashed  line  indicates  the  experimental  value  of  6.55. 
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the  two  finest  grids. 


